
{com}. ****************************************
. **Peter K. Enns and Christopher Wlezien
. **Understanding Equation Balance in Time Series Regression: An Extension
. **Political Science Research Methods
. **Simulation Code to Replicate Table 3
. ****************************************
. 
. *\section[\hspace{c -(}.8in{c )-}Replication Code for Table 3]{c -(}Replication Code for Table 3: Combined Time Series, True Relationship{c )-}
. 
. *********************
. *TABLE 3
. *********************
. 
. *************************
. **T=50; x1, rho=.2, .5, .8
. *************************
. set seed 4545
{txt}
{com}. 
. ******
. *rho=.2
. *****
. *program drop combined_noq
. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 50
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        49    .0250804    .0139806   .0016811    .049245
{txt}
{com}. 
. *****************
. **Column 1 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000    .8939734    .0868375    .406749   1.044851
{txt}{space 7}_b_x1 {c |}{res}      1000    1.001977    .1496983    .566306   1.586026
{txt}{space 6}_sim_3 {c |}{res}      1000   -.8945498    .1709433  -1.473433  -.3026413
{txt}
{com}. 
. *****************
. **Column 2 (%)
. *****************
. *Generate t-statistic for each simulated regression
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000    7.024482    1.462386   2.947782    12.2072
{txt}
{com}. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000    17.18226    8.213022   3.405885   64.90979
{txt}
{com}. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}       998    5.755992    1.474375   2.085295   11.38367
{txt}
{com}. 
. ******
. *\rho=.5
. *****
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 50
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        60    .0260941    .0157092    .000732   .0497449
{txt}
{com}. 
. *****************
. **Column 3 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000    .8946269    .0876222   .4009958   1.063739
{txt}{space 7}_b_x1 {c |}{res}      1000    1.001479    .1402769   .5324845   1.456048
{txt}{space 6}_sim_3 {c |}{res}      1000   -.8954583    .1775381  -1.452201  -.1528786
{txt}
{com}. 
. *****************
. **Column 4 (%)
. *****************
. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000    16.88653    7.770501    3.29589   59.85835
{txt}
{com}. 
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000    6.967871    1.416749   2.651242   11.66704
{txt}
{com}. *
. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}       995    5.731334    1.487475   2.069692   10.73543
{txt}
{com}. 
. ******
. *rho=.8
. *****
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 50
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        63    .0227801    .0131819   .0019436   .0487746
{txt}
{com}. 
. *****************
. **Column 5 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000    .8828905    .0955767   .3714005   1.047931
{txt}{space 7}_b_x1 {c |}{res}      1000    .9968452    .1514654   .6060367   1.574359
{txt}{space 6}_sim_3 {c |}{res}      1000   -.8869969    .1794085  -1.418862  -.3115736
{txt}
{com}. 
. *****************
. **Column 6 (%)
. *****************
. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000    15.96706    7.673462   2.582546   51.83135
{txt}
{com}. 
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000    6.952434    1.505452   2.909203    12.3466
{txt}
{com}. *
. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=2.01

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}       992    5.650685    1.425118   2.059027   10.77631
{txt}
{com}. 
. 
. *************************
. **T=100; x1, rho=.2, .5, .8
. *************************
. 
. ******
. *rho=.2
. *****
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 100
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        58    .0274847    .0138965   .0007935   .0484822
{txt}
{com}. 
. *****************
. **Column 7 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000    .9473803    .0424047   .6869706   1.036117
{txt}{space 7}_b_x1 {c |}{res}      1000    .9991606    .1021567   .7094423   1.292872
{txt}{space 6}_sim_3 {c |}{res}      1000   -.9434747    .1075197  -1.303475  -.5915691
{txt}
{com}. 
. *****************
. **Column 8 (%)
. *****************
. *Generate t-statistic for each simulated regression
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000    9.916342    1.435518   5.217818   14.78982
{txt}
{com}. 
. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000    35.51395    14.58419   9.509529   95.28271
{txt}
{com}. *
. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}      1000    8.933772    1.432699   4.868141   13.51467
{txt}
{com}. 
. ******
. *rho=.5
. *****
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 100
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        53    .0248623    .0140587   .0008648   .0488302
{txt}
{com}. 
. *****************
. **Column 9 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000      .94743    .0440497    .702583   1.032038
{txt}{space 7}_b_x1 {c |}{res}      1000     1.00507    .1039432   .6999887   1.339908
{txt}{space 6}_sim_3 {c |}{res}      1000   -.9520913    .1063337  -1.228035  -.5224339
{txt}
{com}. 
. *****************
. **Column 10 (%)
. *****************
. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000    36.06357    15.74611   9.539879   102.4209
{txt}
{com}. 
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000    9.995309    1.504857   5.648461   15.58384
{txt}
{com}. *
. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}      1000    9.012464    1.453133   5.040073   14.08318
{txt}
{com}. 
. ******
. *rho=.8
. *****
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 100
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        42    .0270249    .0136661   .0056328   .0496919
{txt}
{com}. 
. *****************
. **Column 11 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000    .9428355    .0457272   .6924049   1.025944
{txt}{space 7}_b_x1 {c |}{res}      1000    .9992722    .1059205   .6425325   1.467962
{txt}{space 6}_sim_3 {c |}{res}      1000   -.9436495    .1140591  -1.405189  -.5040581
{txt}
{com}. 
. *****************
. **Column 12 (%)
. *****************
. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000     34.7357    15.10599   9.892268   110.0634
{txt}
{com}. 
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000    9.984527    1.505583   5.829176   15.38713
{txt}
{com}. *
. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}      1000    8.930527     1.45911    3.80729   14.89004
{txt}
{com}. 
. 
. *************************
. **T=5,000; x1, rho=.2, .5, .8
. *************************
. 
. ******
. *rho=.2
. *****
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.2*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. *Simulate the program "combined" N times and save the betas and standard errors.
. *Test whether an equation with mixed orders of integration (combined z, I(0) x1, I(1) x2)
. *can correctly identify TRUE relationships
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        63     .025423    .0143271   .0010209     .04979
{txt}
{com}. 
. *****************
. **Column 13 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000     .998892    .0009148   .9934859   1.000531
{txt}{space 7}_b_x1 {c |}{res}      1000     1.00014    .0147801   .9558638    1.05034
{txt}{space 6}_sim_3 {c |}{res}      1000   -.9991342    .0140588  -1.039089  -.9619758
{txt}
{com}. 
. *****************
. **Column 14 (%)
. *****************
. *Generate t-statistic for each simulated regression
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000    70.71612    1.456268   66.48684   75.33591
{txt}
{com}. 
. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000    1874.251    783.2272   612.7643   5870.345
{txt}
{com}. *
. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}      1000    70.56865    1.414402    65.8921   75.86909
{txt}
{com}. 
. 
. ******
. *rho=.5
. *****
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.5*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        58    .0248132    .0143733   .0004164   .0485559
{txt}
{com}. 
. *****************
. **Column 15 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000    .9989085     .000894    .993852   1.000394
{txt}{space 7}_b_x1 {c |}{res}      1000    .9998737    .0143762   .9562805    1.05159
{txt}{space 6}_sim_3 {c |}{res}      1000   -.9987903    .0141168  -1.043128  -.9596014
{txt}
{com}. 
. *****************
. **Column 16 (%)
. *****************
. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000    1874.488     788.974   676.2466   6100.955
{txt}
{com}. 
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000    70.71539    1.437283   65.50385    75.0206
{txt}
{com}. 
. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}      1000    70.55708    1.451321   65.99915   75.96654
{txt}
{com}. 
. 
. ******
. *rho=.8
. *****
. program drop combined_noq
{txt}
{com}. program define combined_noq, rclass
{txt}  1{com}. drop _all
{txt}  2{com}. set obs 5000
{txt}  3{com}. gen t = _n
{txt}  4{com}. *gen stationary time series (x1)
. gen e1=invnorm(uniform())
{txt}  5{com}. gen x1=e1 if t==1
{txt}  6{com}. replace x1=.8*x1[_n-1] + e1 if t>1
{txt}  7{com}. *gen integrated time series (x2)
. gen u=invnorm(uniform())
{txt}  8{com}. gen x2=u if t==1
{txt}  9{com}. replace x2=x2[_n-1] + u if t>1
{txt} 10{com}. *gen combined times series (z) that his a function of x1 and x2
. gen q=invnorm(uniform())
{txt} 11{com}. gen z = x1 + x2
{txt} 12{com}. tsset t
{txt} 13{com}. reg z l.z x1 l.x1
{txt} 14{com}.         estat bgodfrey
{txt} 15{com}.         mat P = r(p)
{txt} 16{com}.         return scalar pvalue_bg = P[1,1] 
{txt} 17{com}. end
{txt}
{com}. 
. simulate pvalue_bg=r(pvalue_bg) _b _se, reps(1000) nodots: combined_noq
{p2colset 9 19 23 2}{...}

{txt}{p2col :command:}combined_noq{p_end}
{p2colset 1 19 23 2}{...}
{p2col :{txt}[{res:_eq2}]pvalue_bg:}{res:r(pvalue_bg)}{p_end}


{com}. 
. *Percent of simulations which a Breusch-Godfrey test rejects teh null of 
. *no serial correlation
. sum _eq2_pvalue_bg if _eq2_pvalue_bg<0.05

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
_eq2_pvalu~g {c |}{res}        46    .0227756    .0149708   .0002473   .0492189
{txt}
{com}. 
. *****************
. **Column 17 (coef)
. *****************
. sum _sim_1 _b_x1 _sim_3

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 6}_sim_1 {c |}{res}      1000    .9988998    .0009301   .9943075   1.000689
{txt}{space 7}_b_x1 {c |}{res}      1000    1.000551    .0143993   .9526171    1.04721
{txt}{space 6}_sim_3 {c |}{res}      1000   -.9991267    .0145006  -1.038756   -.950473
{txt}
{com}. 
. *****************
. **Column 18 (%)
. *****************
. gen tstat_ly = abs(_sim_1/_sim_5)
{txt}
{com}. sum tstat_ly if tstat_ly>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_ly {c |}{res}      1000     1893.05    809.0085   660.4931   5933.374
{txt}
{com}. 
. gen tstat_x1 = abs(_b_x1/_se_x1)
{txt}
{com}. sum tstat_x1 if tstat_x1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 4}tstat_x1 {c |}{res}      1000     70.7382    1.424628   66.03422   75.09538
{txt}
{com}. 
. gen tstat_lx1 = abs(_sim_3/_sim_7)
{txt}
{com}. sum tstat_lx1 if tstat_lx1>=1.96

{txt}    Variable {c |}       Obs        Mean    Std. Dev.       Min        Max
{hline 13}{c +}{hline 56}
{space 3}tstat_lx1 {c |}{res}      1000    70.55543    1.423106   66.00666   74.93166
{txt}
{com}. 
{txt}end of do-file

